Initial Setup

library(tidyverse)
library(parallel)
library(CoRC)

# helper to run tasks in parallel on all cores
mapInParallel <- function(data, fun, ..., .prep = {}) {
  cl <- makeCluster(detectCores())
  clusterCall(cl = cl, fun = eval, .prep, env = .GlobalEnv)
  ret <- clusterApplyLB(cl = cl, x = imap(data, ~ set_names(list(.x), .y)), fun = lapply, as_mapper(fun), ...)
  stopCluster(cl)
  flatten(ret)
}

3D Trajectory Plot of a Calcium Model

This example loads the Kummer2000 - Oscillations in Calcium Signalling model. The model has 3 species which oscillate. These oscialltions can be visualized as a trajectory through a 3D space. The example does this once in a deterministic and once in a stochatic fashion.

loadSBML("http://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000329")
#> # A COPASI model reference:
#> Model name: "Kummer2000 - Oscillations in Calcium Signalling"
#> Number of compartments: 1
#> Number of species: 3
#> Number of reactions: 8

# Run 24 sec (2 Periods)
setTimeCourseSettings(24, intervals = 10000)

timeseries <- list(
  deterministic = runTimeCourse()$result,
  stochastic = runTimeCourse(method = list(method = "directMethod", use_random_seed = T, random_seed = 1))$result
)

# simplify species names so they are valid R syntax
timeseries <- map(timeseries, set_tidy_names, TRUE)

# Create the same plot for both timeseries
plots <- map(
  timeseries,
  plotly::plot_ly,
  type = "scatter3d",
  mode = "lines",
  x = ~ G.alpha,
  y = ~ activePLC,
  z = ~ Calcium,
  color = ~ Time
)

unloadModel()

plots$deterministic

plots$stochastic

Statistics of Repeated Stochastic Simulations

This implements an example from the Condor-COPASI paper. The example illustrates advantages of parallel processing.

# Run 1000 stochastic time series possibly in parallel
loadModel("https://raw.githubusercontent.com/copasi/condor-copasi/master/examples/stochastic_test.cps")
#> # A COPASI model reference:
#> Model name: "Kummer calcium model"
#> Number of compartments: 1
#> Number of species: 3
#> Number of reactions: 8

# timeseries <- 1:1000 %>% map(~ runTimeCourse()$result)
timeseries <-
  # Defines parallel evaluation:
  mapInParallel(
    # perpare worker (.prep),
    .prep = quote({
      library(CoRC)
      loadModel("https://raw.githubusercontent.com/copasi/condor-copasi/master/examples/stochastic_test.cps")
      setTimeCourseSettings(method = list(method = "directMethod", use_random_seed = T))
    }),
    # iteration data (1000 random seeds),
    1:1000,
    # iteration code (formula ~)
    ~ runTimeCourse(method = list(random_seed = .x))$result
  )

# Combine all results and reshape the data
plotdata <-
  timeseries %>%
  bind_rows() %>%
  group_by(Time) %>%
  # calculate mean and sd for all time points
  summarise_all(funs(mean, sd)) %>%
  # gather all values so the column "name" identifies "a_mean", "b_sd" etc.
  gather("name", "value", -Time) %>%
  # split up information on species (a,b,c) and type of value (mean, sd)
  separate(name, c("species", "type"), "_") %>%
  spread(type, value)

print(plotdata, n = 6)
#> # A tibble: 2,403 x 4
#>     Time species  mean    sd
#>    <dbl> <chr>   <dbl> <dbl>
#> 1 0.     a        8.00 0.   
#> 2 0.     b        8.00 0.   
#> 3 0.     c        8.00 0.   
#> 4 0.0500 a        7.06 0.249
#> 5 0.0500 b        8.12 0.117
#> 6 0.0500 c        5.60 0.442
#> # ... with 2,397 more rows

plot <-
  ggplot(data = plotdata, aes(x = Time, y = mean, group = species, tt_sd = sd)) +
  geom_ribbon(aes(ymin = mean - sd, ymax = mean + sd, fill = species), alpha = 1 / 4) +
  geom_line(aes(color = species)) +
  guides(fill = "none") +
  labs(
    x = paste0("Time (", getTimeUnit(), ")"),
    y = paste0("Concentration (", getQuantityUnit(), ")"),
    color = "Species"
  )

unloadModel()

plotly::ggplotly(plot, tooltip = c("group", "x", "y", "tt_sd"))

Parameter Scan

2D Scan Over the Cartesian Product of Two Species Concentration Vectors

This implements an example from the Mendes2009 paper on COPASI use cases.

loadSBML("https://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000068")
#> # A COPASI model reference:
#> Model name: "Curien2003_MetThr_synthesis"
#> Number of compartments: 1
#> Number of species: 9
#> Number of reactions: 3

setSteadyStateSettings(calculate_jacobian = FALSE, perform_stability_analysis = FALSE)

# Cartesian product of the input values
scan <- cross_df(
  list(
    cysteine = 0.3 * 10 ^ seq(0, 3, length.out = 6),
    adomed = seq(0, 100, length.out = 51)
  )
)

scan <-
  scan %>%
  mutate(
    # Calculate steady state fluxes for every row
    ss_fluxes = pmap(., function(cysteine, adomed) {
      setSpecies("Cysteine", initial_concentration = cysteine)
      setSpecies("adenosyl", initial_concentration = adomed)
      ss <- runSteadyState()
      stopifnot(ss$result == "found")
      ss$reactions$flux
    }),
    # The second flux is CGS
    CGS = map_dbl(ss_fluxes, 2),
    # The third flux is TS
    TS = map_dbl(ss_fluxes, 3)
  )

plot <-
  ggplot(data = scan, aes(x = adomed, group = cysteine)) +
  geom_point(aes(y = CGS, color = "CGS")) +
  geom_point(aes(y = TS, color = "TS")) +
  geom_line(aes(y = CGS, color = "CGS")) +
  geom_line(aes(y = TS, color = "TS")) +
  labs(
    x = paste0("Adomed (", getQuantityUnit(), ")"),
    y = paste0("Flux (", getQuantityUnit(), " / ", getTimeUnit(), ")"),
    color = "Species"
  )

unloadModel()

plotly::ggplotly(plot)

2D Scan over Random Concentrations of Two Species

This implements an example from the Mendes2009 paper on COPASI use cases. It is in many ways similar to the previous example but is written to run parallelized.

loadSBML("https://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000068")
#> # A COPASI model reference:
#> Model name: "Curien2003_MetThr_synthesis"
#> Number of compartments: 1
#> Number of species: 9
#> Number of reactions: 3

# 10000 repeats of steady state task with random cysteine and adomed
scan <-
  # Defines parallel evaluation:
  mapInParallel(
    # perpare worker (.prep),
    .prep = quote({
      library(CoRC)
      loadSBML("https://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000068")
      setSteadyStateSettings(calculate_jacobian = FALSE, perform_stability_analysis = FALSE)
    }),
    # iteration data (10000 random seeds),
    1:10000,
    # iteration code (formula ~)
    ~ {
      set.seed(.x)
      cysteine <- 0.3 * 10 ^ runif(1L, min = 0, max = 3)
      adomed <- runif(1L, min = 0, max = 100)
      setSpecies(
        key = c("Cysteine", "adenosyl"),
        initial_concentration = c(cysteine, adomed)
      )
      ss <- runSteadyState()
      stopifnot(ss$result == "found")
      list(
        cysteine = cysteine,
        adomed = adomed,
        CGS = ss$reactions$flux[2],
        TS = ss$reactions$flux[3]
      )
    }
  )

# Combine all results and reshape the data
plotdata <-
  scan %>%
  bind_rows() %>%
  gather("reaction", "flux", CGS, TS)

plot <-
  ggplot(data = plotdata, aes(x = adomed, y = flux, group = reaction, tt_cys = cysteine)) +
  geom_point(aes(color = reaction), alpha = 1 / 10, size = 3 / 4) +
  labs(
    x = paste0("Adomed (", getQuantityUnit(), ")"),
    y = paste0("Flux (", getQuantityUnit(), " / ", getTimeUnit(), ")"),
    color = "Species"
  )

unloadModel()

plotly::ggplotly(plot, tooltip = c("tt_cys", "x", "y"))

Parameter Estimation